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The influence of uncorrelated, quenched disorder on the phase transition of 
two dimensional Potts models will be reviewed. After an introduction where 
the conditions of relevance of quenched randomness on phase transitions are 
exemplified by experimental measurements, the results of perturbative and 
numerical investigations in the case of the Potts model will be presented. 
The Potts model is of particular interest, since it can have in the pure case 
a second-order or a first-order transition, depending on the number of states 
per spin. In 2D, transfer matrix calculations and Monte Carlo simulations 
were used in order to check the validity of conformal invariance methods 
in disordered systems. These techniques were then used to investigate the 
universality class of the disordered Potts model, in both regimes of the pure 
model phase transitions. A test of replica symmetry became possible through 
a study of multiscaling properties and a detailed analysis of the probability 
distribution of the correlation functions was also made possible. 



1 Introduction 



Quenched disorder has been the subject of an intensive activity in statistical physics 
during the past decades. The quahtative influence of disorder coupled to the energy- 
density at second order phase transitions is well understood since Harris proposed a cel- 
ebrated relevance criterion At first order transitions, randomness obviously softens 
the transitions, and, under some circumstances may even induce second order transitions 
according to a picture first proposed by Imry and Wortis Q and then stated on more 
rigorous grounds by Aizenman and Wehr ||3|, Q , implying in particular that an infinites- 
imal disorder induces continuous transitions in 2D. Reviews can be found in the books 
of S.K. Ma, J.L. Cardy or Vik. Dotsenko |, |, Q. 

In spin models, the influence of quenched disorder strongly depends on the nature 
of randomness, i.e. to which quantity the perturbation is coupled in the Hamiltonian. 
Consider for example a very general model with spins s, on a lattice and define the 
Hamiltonian 

-pn = J2 KijSiSj + J2 HiS, + J2 DisiUif + ... 

(ij) i i 

where i^y, Hj, or rij are independent random quenched variables drawn from some prob- 
ability distributions P[Kij\, P[Hj], or -P[nj], and which respectively describe random- 
bond or dilute problems ||8|, random fields |^, 0, and random anisotropy 



models |T5|. As usually, the sum over (ij) is supposed to be restricted to nearest 



neighbours. Usually, we decide to work with uncorrelated quenched random variables, 
for example Kij = J KV[K]dK = Kq, and KijKf^i = AdikS^i and we will here only con- 
centrate on the first category, namely random-bond systems where disorder is coupled to 
the energy density. Special cases of probability distributions of interest are for example 
(all written here in the bond version of the problem) 

i) dilution problems, where non magnetic impurities are randomly distributed on the 

bonds or sites of the lattice, e.g. in the bond case 

nK^J] = UlP^iK^j -K) + {l-p)5{K,j)l (1) 

ii) binary distributions, where we can imagine for example a disordered alloy of two 

magnetic species 

nKi,] = X{\p5{Ki, -K) + {1- p)6{Kij - Kr)], (2) 

iii) Gaussian distributions which are of particular interest to perform Gaussian integra- 

tion in analytic approaches, 

'Pi^i^] = Ili-7=^ exp (-(K,, - Kf/2a'')] (3) 
{ij) ^^'^^ 
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iv) Continuous self-dual distributions (to be discussed later) 



K,, 



(4) 



In each case, there is a control parameter (p, r, a or A) which determines the strength of 
disorder. The phase diagram is sketched in figure [l| for dilution and binary distribution. 
We expect a transition line between ordered and disordered phases along which the 
transition is continuous in 2D. For dilute problems, there exists a percolation threshold 
where the transition temperature vanishes and below which long range order cannot 
exist. In the bimodal case, the percolation fixed point at finite q is reached in the limit 
r — > oo. On the transition line, there should be some particular strength of disorder 
corresponding to the location of the random fixed point, where corrections to scaling 
should be small. 
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Figure 1: Phase diagram of the dilution (left) and binary problems (right). 



The g— state Potts model [16| is the natural candidate for the investigations of in- 
fluence of disorder, since the pure model exhibits two different regimes (see figure §): a 
second order phase transition when g < 4 and a first order one for g > 4 in two dimen- 
sions {2D). In 3D, ordering is easier and the transition becomes weakly first-order at 
(? = 3 already. 

The 2-dimensional g-state Potts model is defined by the following Hamiltonian : 



(5) 



where the sum is restricted to nearest neighbours (here on a square lattice), the de- 
grees of freedom {ai} can take q values and the exchange couplings Kij = Jij/ksT are 
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Figure 2: Phase diagram of the 2D pure Potts model as a function of the number of states 
per spin. The transition is of second-order when g < 4 and becomes of first- 
order above. Inserts show Monte Carlo snapshots of typical spin configurations. 



quenched independent random variables chosen according to some distribution {Kij} to 
be specified later. 

Many results were obtained for quenched randomness in this model in the last ten 
years, including 



) perturbative expansions for 2 < g < 4 |17|, |]J, |19|, |20|, |2g, g3|, |2j, g5||26|, g7| . This 



work was initiated in a series of papers by A.W.W. Ludwig and J.L. Cardy [17, 18 
19]. A systematic study of energy-energy and spin-spin correlators in random bond 



Ising and Potts models, including tests of replica symmetry, was then performed 
by a group around Vl.S. Dotsenko pOl, |2l|, 0, 0, |4|, EBI, with M. Picco, P. 



Pujol, Vik. Dotsenko, J.L. Jacobsen, and M.-A. Lewis. 

ii) Monte Carlo simulations in both regimes 2 < q < A and q > A 

32, 33, 34, |3^. In the second order regime, simulations were first performed by 
S. Wiseman and E. Domany in the case of the Ashkin- Teller model |3^] (which 
exhibits a critical line and for a given set of couplings belongs to the 4— state Potts 
model universality class) and in the 3— state Potts model by J.K. Kim [37|. In 
the case q = 8, the first Monte Carlo simulations were reported by S. Chen, A.M. 
Ferrenberg and D.P. Landau in Refs. 39 1, but then in both regimes, more 
accurate results were obtained by different groups, like M. Picco 



T. Olson 



and P. Young [34|, or the present authors |29 



iii) Transfer matrices in both regimes [^, 41, |4^, 43, 44, 46, 47, 48 1. This technique 
was used to study random Potts models very early by U. Glaus, then more refined 



4 



computations were reported, for example by J.L. Cardy and J.L. Jacobsen 43 1 
or in Ref. 



, were 



iv) High-temperature series expansions, initially used in the random Ising model [49, 50 1 

to show the logarithmic corrections, then extended at g = 3 in 3D |51, 52 
also shortly applied to the two-dimensional case ||^ 

v) Short-time dynamic scaling |54, 56, 57, |5^ . 

We also notice that dynamical properties (in order to discriminate between conventional 
or activated dynamics) linked to non self-averaging have been recently studied also [59| 
and that the interesting limit g — > oo was carefully investigated in Refs. [45, Reviews 
on selective parts of the subject were reported in Ph.D dissertations in Refs. [61, |65, |63 



Although closely related, the random-bond Ising model will not be discussed here. 

We will only remind 



since it was already the subject of many reviews, e.g. [66, 67, |65 
here that in the pure 2D Ising model, the exponent a of the specific heat vanishes and 
further investigation is needed, since the Harris criterion becomes inconclusive. It is 
now generally believed that uncorrelated quenched randomness is a marginally irrele- 
vant perturbation which does not modify the universal critical behaviour, but produces 
logarithmic corrections. For example the singular part of the specific heat exhibits the 
following behaviour. 



Cs{t)^Hi/\t\) 

C,(t)~lnln(l/[t[) 



9<.\t\<. 1, 

\t\ « e, 



~ e 



where g is the amplitude of disorder, linked to the impurity concentration. The behaviour 
of the main physical quantities in the neighbourhood of the random fixed point is given 
in table |^. It was confirmed unambiguously by Monte Carlo simulations [69, 70, 71, 72 [ 
as shown in the same table. We also note that conformal mappings associated to Monte 
Carlo simulations were initially used in the random-bond Ising problem by A. Talapov 



and Vl.S. Dotsenko [73[ 



Random-bond Ising model 



analytical results numerical results 



Correlation function 






1 = 


0.2493 ±0.0014 


Correlation length 




\t\-\Hi/\t\)Y/^ 






Magnetisation 


m{t) ^ 


\t\'^\Hi/\t\)Y"' 


P/v 


= 0.1245 ±0.0009 


Susceptibility 


x{t)- 






= 1.7507 ±0.0014 


Specific heat 




lnln{l/[t[) 







Table 1: Critical behaviour of the random-bond Ising model. 
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The numerical studies of disordered models showed that many difficulties, like the 
lack of self averaging |74, 75, 76, ^ or varying effective exponents due to crossover 



phenomena. Averaging physical quantities over the samples with a poor statistics may 
thus lead to erroneous determinations of the critical exponents. Almost all the studies 
mentioned here were reported in the case of the random bond system with self-dual 
probability distributions of the coupling strengths in order to preserve the exact knowl- 
edge of the transition line, which is an important simplification when one wants to use 
finite-size-scaling techniques or conformal mappings which hold at criticality only. 

In real experiments on the other hand, disorder is inherent to the working-out process 
and may result e.g. from the presence of impurities or vacancies. For the description 
of such a disordered system, dilution is thus more realistic than for example a random 
distribution of non- vanishing couplings (the so-called random-bond problem). Since 
universality is expected to hold, the detailed structure of the Hamiltonian should not 
play any determining role in universal quantities like critical exponents, but crossover 
phenomena may alter the determination of the universality class. Experimentally, the 
role of disorder in 2D systems has been investigated in several systems. Illustrating 
the infiuence of random defects in the case of the 2D Ising model universality class, 
samples made of thin magnetic amorphous layers of (Tbo.27Dyo.73)o.32Feo.68 of 10 A 
width, separated by non magnetic spacers of 100 A Nb in order to decouple the magnetic 
layers were produced using sputtering techniques. A structural analysis (high resolution 
transmission electron microscopy and x-ray analysis) was performed to characterise the 
defects inherent to such amorphous structures, and in spite of these random defects 
separated on average by a distance of a few nm, the samples were shown to exhibit 
Ising-like singularities with critical exponents 



/? = 0.126(20), 7 = 1.75(3), (^=15.1(10). 

This is coherent with the fact that disorder does not change the universality class of the 
2D Ising model, apart from logarithmic corrections which are probably impossible to 
observe experimentally, since their role becomes prominent only in the very neighbour- 
hood of the critical point. A beautiful experimental confirmation of the Harris criterion 
- which predicts a modification of the critical behaviour in random systems when the 
exponent a of the specific heat is positive for the pure system - was reported in a Low 
Energy Electron Diffraction investigation of a 2D order disorder transition ||7^ belong- 
ing to the 4-state Potts model universality class. Order-disorder transitions of adsorbed 
atomic layers are known to belong to different two-dimensional universality classes de- 
pending on the type of superstructures in the ordered phase of the ad-layer |8^, |8^. The 
substrate plays a major role in ad-atom ordering, as well as the coverage (defined as the 
number of ad-atoms per surface atom) which determines the possible superstructures of 
the over-layer. For example, sulfur chemisorbed on Ru(OOl) exhibits four-state or three- 
state Potts critical singularities for the p{2 x 2) and the (\/3 x \/3)i?30° respectively 
(at coverages 1/4 and 1/2). The case of the (2 x 2)— 2H/Ni(lll) order-disorder transi- 
tion of hydrogen adsorbed on the (111) surface of Ni thus belongs to the 2D four-state 
Potts model universality class, since the ground state, stable at low temperatures, has 
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a four-fold degenaracy due to the four possible coverings of the ad- atoms at the (111) 
surface (see figure ^) . 




Figure 3: (2 x 2)— 2H/Ni(lll) order-disorder transition of hydrogen. The ground-state 
has a four-fold degenaracy due to the four possible covering of the (111) surface 
of Ni by H ad-atoms. 




Figure 4: (2 x 2)— 2H/Ni(lll) order-disorder transition of hydrogen with oxygen impu- 
rities randomly chemisorbed on the surface, and occupying some of the lattice 
sites. 

The expected exponents are thus close to theoretical values of (3 = 1/12 ~ 0.083, 
7 = 7/6 ~ 1.167, and = 2/3 ~ 0.667 for example. 

Low energy electron diffraction (LEED) makes possible to measure these exponents 
through the diffracted intensity /(q) or structure factor. This is the two-dimensional 
Fourier transform of the pair correlation function of ad-atom density. Long range fluc- 
tuations produce an isotropic Lorentzian centered at the superstructure spot position qo 
with a peak intensity given by the susceptibility and a width determined by the inverse 
correlation length, while long range order gives a background signal proportional to the 
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order parameter squared: 



I(q) = (m^)J(q- qo) + ^ 



l + e2(q-qo)2- 
The following exponents were thus measured |8^, Q 

fi = 0.11 ± 0.01, 7 = 1.2 ± 0.1, u = 0.68 ± 0.05 

in correct agreement with 4-state Potts values (the small deviation, especially for the ex- 
ponent fj, is attributed to the logarithmic corrections to scaling of the pure 4-state Potts 
model [ p^ ] ) . The same experiments were then reproduced in the presence of intentionally 
added oxygen impurities, at a temperature which is above the ordering temperature of 
pure oxygen adsorbed on the same substrate. The mobility of these oxygen atoms is fur- 
thermore considered to be low enough at the hydrogen order-disorder transition critical 
temperature that they essentially represent quenched impurities randomly distributed 
in the hydrogen layer (see figure Q). The exponents become 

(3 = 0.135 ± 0.010, 7 = 1.68 ± 0.15, v = 1.03 ± 0.08, 

which definitely supports the modification of the universality class in the presence of 
quenched disorder, in agreement with Harris criterion {a = 2/3 for the 4-state Potts 
model). 

The aim of this lecture is to provide a review of both perturbative and numerical 
studies of the bond diluted Potts model for several values of the number of states per spin 
(in order to cover the two different regimes of the pure system's phase transitions). In 
section ^, the perturbative approach is discussed and the essential results are summarised 
in the 2D case with 2 < g < 4. The details of the numerical techniques used in two 
dimensions are presented in section ^ and |^, while the comparison between numerical 
and analytical results is made in |^. 

2 Perturbative approach in 2D 
2.1 Replicas and relevance criterion 

For a specific disorder realization [Kij], the Hamiltonian T-L[Kij,ai\ is written 

- m\Kij,o,] = J2(Ko + SK,j)6^^,^^ . (6) 

The corresponding partition function and free energy are given as follows 

Z[K^j] = J P[ai]e-^^[^'^''^'l, 
F[Kij] = -kBT In Z[K,j]. 
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To get the quantities of interest, for example the average free energy, we have to perform 
an average over the distribution V[Kij]^ 

F = F[K~] = -ksT J V[Kij]V[K^j] In Z[K^j]. (7) 

Averaging the logarithm of the partition function is possible through the identity 

lnZ= lim 

n—>oo n 

which requires, before averaging, to produce n copies (with labels a) of the system with 
the same disorder configuration, 

and then to perform integrations over VlKij], 



-X - ^-x+^{x^-x^)+.. 



Q-yy — g " i 2 



leading to 



(zmr = / n ^[^1 



\a=l 



xe I j ^ 3 



In the leading term, 6Kij has RG eigenvalue yt = d—Xg, and corresponds to a simple shift 
of the transition 4 temperature (which is obviously a relevant effect). In the next term, 
the second moment of the distribution, SK"^ — 5K , has RG eigenvalue yn = d — 2x£, and 
all the following terms are irrelevant^. Using hyperscaling relation, the Harris scaling 
dimension of disorder is rewritten 

yjj = a/v. (9) 

It implies that at second-order transitions, disorder is a relevant perturbation which 
modifies the critical behaviour when the specific heat exponent a of the pure system 



^In the case of an annealed disorder, the impurities are thermahzed (this would only be possible 
if the relaxation time of randomness is small compared to the time scale of the experiment) and 
their probability distribution VlKij] depends strongly on the spins (and vice-versa); since it is the 
equilibrium distribution V[Kij] = J T)[ai]Z~^e~'^''^^^'-^''''\ where Z does no longer depend on the 
disorder realization, but is obtained through Z = J 'D[Kij]Ty[ai]e~^'^^'^'^ ''^^K In the annealed case, if 
the impurity concentration is kept constant, there is a 'Fisher's renormalization' of the exponents if 
the specific heat of the pure system is diverging p^ . 

^The leading (unperturbed) term is written in the continuum limit as — /37ic ~ rno f £a(r)d^r 
where mo stands for Ko + 5K while the perturbation is written go J '^a^p ^a{T)ei3{r)d^r with go 
corresponding to 5K^ — SK^ . 
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is positive, while it is irrelevant (and universal properties are thus unaffected by ran- 
domness) when a is negative. In the borderline case a = 0, randomness is marginal 
to leading order. This is the case for example of the ID Ising model discussed in the 
introduction, where quenched disorder is eventually marginally irrelevant and produces 
only logarithmic corrections to the unchanged leading critical behaviour. 

The case of first-order transitions was considered later, by Imry and Wortis, Aizen- 
man and Wehr, then Hui and Berker g |, |]. It can be intuitively understood from 
the above results simply by noticing that the existence of a latent heat at first-order 
transition corresponds to a discontinuity of the energy density and can be described by 
a vanishing energy density scaling dimension, so that disorder is always relevant in this 
sense. 



2.2 Perturbation techniques 
2.2.1 Average correlation functions 

Many results were obtained in the 2D random Potts model using RG perturbations, 
mainly around Ludwig and VI. Dotsenko. In equation (|8|), it appears that couples 

the replicas via energy-energy interactions Yliai^iii.^^'^ ~ ) Z^r £a(r)e/3(r) which have 
to be treated as a perturbation around the pure fixed point. Here, eQ.(r) is a short 
notation for 5 (a) and r stands for the lattice unit vectors. The second cumulant 

a- a ■ 

» 3 

of the coupling distribution will be denoted (70 in the following. Two different schemes 
were considered in the literature |25, Sl|, 



i) replica symmetric scenario, where all the replicas are coupled through the same in- 

teraction strength, 

5^ 50 E^- We/5 (10) 

ii) replica symmetry breaking scenario, where the coupling between replicas are replica- 

dependent, 

E 5a/3Eea(r)e/3(r). (11) 

The program is thus to consider 2D Potts model with weak bond randomness, compute 
the scaling dimensions x'^{n) of the order parameter and x'^{n) of the energy-density 
perturbatively^ around the Ising model conformal field theory, and then take the replica 
limit n ^ 0. Expansions are performed around the pure model fixed point (weak dis- 
order) in terms of the disorder strength ^K^- — 6Kij , and the exponents are given in 
powers of i/h = ajv. 

Different expansion parameters can be found in the literature, and it is worth col- 
lecting the main notations. The Potts models can be identified to minimal conformal 

^The primes denote the scaling dimensions at the random fixed point. 
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models [87] which are parametrised by an integer m which determines the central charge 
and critical behaviour of the model. The correspondence is given by 

vr 

= n — ^ , N — 1 

cos~i(^/g/2) 

and the central charge and exponents follow from 

c = 1 



m{m + 1) ' 
m + 3 
^ 2m ' 

(m + 3)(m — 1) 
8m{m + 1) 

We note that m = 3 for the Ising model, m = 5 for the 3— state Potts model and m — > oo 
for the 4— state Potts model. The Harris RG eigenvalue becomes 

m — 3 



m 

which is proportional to g — 2 to linear order {q being the number of states per spin of 
the Potts model): 

yj, = ±{q-2)-^{q-2f + 0[{q - 2)% (12) 

Using a Coulomb gas representation, a natural expansion parameter e is defined through 

2 m + 1 4 
mo 

and it is linked to yn hy e = — ^yn- 

i) The replica symmetric case is based on the assumption that replica symmetry is not 
broken initially, and is then preserved by the renormalization group. The renor- 
malization of the coupling constant qq is determined by perturbative calculation 
using the operator algebra. For any scaling operator (j), the perturbed two-point 
correlation function {(p{0)(f)(R))g corresponds, in the limit n — > 0, to the average 
correlator {(p(0)(j)(li)) . We can write 

Tr0(O)(^(R)e-^(^=+^'') 



where the perturbation term —fSTig = go J J2a^i3 £a(r)£/3(r)d^r acts on the 'critical' 
Hamiltonian —f3Ti.c = ruQ J J2a^a{^)d'^r + ho JJ2a'^ai^)d'^''^^ ^^^t term being 
included in order to compute the renormalization of the spin operator. 
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When expanded in terms of unperturbed correlators, it yields the following expan- 
sion 1881, |89|1 



1 

2' 



The renormalization of the coupling constant follows from the expansion [20, 21 



23, ill 



90 ea(r)e/3(r)dV 

+\gljY.^-{r)ep{v)d\ J Ye,{r')es{r')d\' + ... 



= g j Ye^{r)ep{v)d\ (13) 

leading to 5 = (70(1 + -^i9o + ^2^0 + •••)• The successive terms are obtained from 
operator product expansions. In the first order term, the dominant contribution to 
Ai follows from contraction of neighbouring pairs, e(r)e(r') ~ |r — r'|~^^=, in the 
same replica (/? = 7 7^ a, 5 and a ^ 5). Such an expression has to be understood 
inside unperturbed correlators. Including combinatorial factors (they are 2(n — 2) 
such factors), integration over space up to an infrared cutoff, h> |r — r'|, leads to 
^OqAi / £o(r)£5(r)d^r, where Ai is dominated by 

Ai = 2{n - 2) / |r - r'|-2^-dV' = 47r(n - 2)- 



|r-r'|<fe 2 - 2x, 



£ 



Since yn = 2 — 2xe, one recovers the Harris criterion. Up to the first order, 
we get the following expression for g in terms of the bare coupling constant go, 
g = 5(0(1 + 27r(n — 2)-^by^gQ + OigK))- Following Dotsenko and co-workers, the 
coupling constants are multiplied by a factor in order to get dimensionless 
coupling constants g{h). 

The /?— function up to second order in the n — > limit is finally given by 

m = ^ = yngib) - Sng^b) + 327rV(6) + OigHb)). (14) 
It leads to a non trivial IR fixed point (which determines the long distance physics) 

5c = ^y/f + ig^y|^ + 0(yf^). 

Renormalization of the energy and the order parameter density operators follow 
from the same analysis, e.g.: 



mo / E^-(r)dV (1 + 50 / E e/3(r')^7(r')dV 
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+\9lJ E^/3(r'K(r')dV'| ^e,(r")e,(r")dV + ... 
= m ^ea(r)d^r, (15) 



and provide the expansions (details of the calculation of the integrals, using a 
Coulomb gas representation, can be found e.g. in Ref. |^0[) m = mo(l + Bigo + 
-B250 + • ■ ■) = ^eW-o and h = ho{l + Cigo + (^250 + • • •) = ^o-Zio, leading when 
n — > to 

ls{g) = ^ = -47.5(&)+8vrV(6), 

7.(.) ^ ^ - W(5) (1 + Vv'I'iS) ) ^ 

(16) 

For the correlators themselves, the renormalization equations can be written 

(<^(o)<^(.R)) = -^^^s-^'^^imm)), 

where s is the scaling factor. Using now 



dlnZ,^ 
= ^hTT' 

or InZ^ = / 'j^{g)d\nb, the ratio z\g{h)) rewritten 

which is dominated at long distances by g ~ gc, such that /^'''^7(^(g)dln6 ~ 
70(g'c)liis. The homogeneity equation thus becomes 

((/.(O)0(sR)) = s-2(^^-T^(sc))^^(o)0(R)), 

and choosing a rescaling factor s = , the two point correlator decays as 

((/)(O)0(R)) ~i?-2{x^^749c))_ (17) 

The corresponding scaling dimension is modified according to 

x'^ = x^--i4>{gc)- (18) 
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Collecting the results of Dotsenko and co-workers, we give the new thermal and 
magnetic scaling dimensions (with primes) in terms of the original ones (un- 
primed) [ 



21, 22, 2 



Xe - leig, 
1 



Xa + 



1 r^(-§)r2(i) 



32r2(-i)r2(-i 



(19) 
(20) 



ii) The other assumption of a broken replica symmetry leads to a different fixed point 
structure. The coupling between replicas, gap, is now dependent on the pair in- 
dexes, and it is generalised to a continuous variable x instead of pair indices, 
ga/3 ~^ 9{oi — (3) = g{x). It is found that there is only one marginally attractive so- 
lution for the coupling g{x) which then enable to compute 7^(5') and ja{g), leading 
to a modified thermal exponent 

Xe=Xs + lyH + 0{yH), 



21 



while to yjj order, the magnetic scaling index remains the same as in the replica 



symmetric scenario |22]. 



2.2.2 Multiscaling and higher order moments of the correlators 

In order to measure some other differences between the replica symmetric and the replica 
symmetry breaking cases, the moments of the correlators can also be helpful. For any 
scaling field </>(r), multiscaling arises when the scaling dimensions associated to the 
moments of the correlators do not follow a simple linear law: 

(0(O)^(R))?^i?-2px^'', 



(22) 



In the magnetic sector, a difference between the two cases occurs to order for the 
second moment |25]: 



RS =x^-lyH + ^(41n2-ii)y2^ + 0(y|,), (23) 

RSB x';. =x.-ly^^ + l(41n2-iL)y|, + o(y3,). (24) 

For higher order moments, the computation were only performed in the replica symmetric 
case 1 19, 26, 27], leading to 



97r2 
1 



(3p-4)iq-2f+0[{q-2Y 



16 



32 



iP - l)[i|-41n2+i(33-29v^V3){p-2)]y|^ + 0{yfl) 



(25) 



(26) 



14 



2.2.3 Are these effects measurable? 

The question is now to try to detect numerically the effects discussed above. These are 
perturbative expansions around q = 2, so that a natural choice of system to measure 
the scaling dimensions in the presence of quenched disorder is the 3— state Potts model. 
At (; = 3 we have = 4/5 and yn = 2/5 from which the perturbed scaling dimensions 
in the energy and magnetic sectors can be obtained. The values are given in table g 
The numerical data clearly show that the expected variations are quite small and need 
accurate numerical techniques to discriminate between RS and RSB schemes. 



Scheme 




Scaling dimensions 
















Pure system 


0.13333 


0.800 


0.13333 


0.13333 


0.800 


RS 


0.13465 


1.000 


0.11761 


0.18303 


1.090 


RSB 


0.13465 


1.020 


0.12011 







Table 2: Comparison between pure 3— state Potts model critical exponents and the ex- 
pected values obtained from perturbation expansions. The notation x^2 corre- 
sponds to the second moment of the spin-spin correlation function, while x^o 
and are associated to the typical correlations. 
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3 Numerical techniques in 2D 



3.1 Monte Carlo simulations 
3.1.1 Cluster algorithms 

For the simulation of spin systems, standard Metropolis algorithms based on local up- 
dates of single spins suffer from the well known critical slowing down. As the second-order 
phase transition is approached, the correlation length becomes larger and the system con- 
tains larger and larger clusters in which all the spins are in the same state. Statistically 
independent configurations can be obtained by local iteration rules only after a long 
dynamical evolution which needs a huge number of MC steps. This makes this type of 
algorithm inefficient close to a critical point. 

Since the transition of the disordered Potts model is always expected to be on average 



a second-order one, the resort to cluster update algorithms is more convenient |90, 91 1. 
The main recipe of cluster algorithms is the identification of clusters of sites using a bond 
percolation process connected to the spin configuration. The spins of the clusters are 
then independently flipped. A cluster algorithm is particularly efficient if the percolation 
threshold coincides with the transition point of the spin model, which guarantees that 
clusters of all sizes will be updated in a single MC sweep. 

In the case of the Potts model, the percolation process involved is obtained through 
the mapping onto the random graph model. These algorithms are based on the Fortuin- 
Kasteleyn representation [B2| where bond variables are introduced. In the Swendsen- 



Wang algorithm ||9^, a cluster update sweep consists of three steps: depending on 
the nearest neighbour exchange interactions, assign values to the bond variables, then 
identify clusters of spins connected by active bonds, and eventually assign a random 



value to all the spins in a given cluster. The Wolff algorithm |94] is a simpler variant 
in which only a single cluster is flipped at a time. A spin is randomly chosen, then 
the cluster connected with this spin is constructed and all the spins in the cluster are 
updated. 

Both algorithms considerably improve the efficiency close to the critical point and 
their performances are comparable in two dimensions, so in principle one can equally 
choose either one of them. Nevertheless, when one uses particular boundary conditions, 
with fixed spins along some surface for example, the Wolff algorithm is less efficient, 
since close to criticality the unique cluster will often reach the boundary and no update 
is made in this case. 



3.1.2 Definition of the physical quantities 

For each disorder strength, many samples (A'^rdm) from the same probability distribution 
are studied at a given temperature. Each sample, initialised from the low-temperature 
phase, is thermalized during A'th Monte Carlo sweeps and the physical quantities are 
then measured during A'^mc sweeps. Different quantities, averaged over the MC sweeps 
denoted by (...), are conserved for all the samples. Many physical quantities can be 
measured: 
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i) The order parameter density follows from the standard definition for the Potts 
model: 

Ml \ QPmax 1 

q-1 

where Pmax is the fraction of spins in the majority orientation 

Pmax = MaXa (pa), Pa = ^ ^^a^ ,a • 

j 

Thermal average is understood in the notation M. To obtain the local order 
parameter (cr(i)) at site i, it is counted 1 when the spin at site i is in the majority 
state and otherwise. 

ii) The susceptibility is given by fluctuation-dissipation theorem 

kBTx = {a^)-{af. 

iii) Energy density: 

iv) Specific heat: 



C/kB = - {ef. 
v) Energy Binder cumulant: 

iv) Correlation functions: the connected spin spin correlation function Ga{i-ij) = 
{(y{i)cr{j)) — (f^) at criticality is obtained by the estimator of the paramagnetic 
phase, 

- 1 

given by the probability that spins at sites i and j belong to the same finite cluster. 
All these quantities are then averaged over the disorder realisations 

V)= j{...)v[{..m-)- 
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3.2 Transfer matrix technique 



The disordered Potts model can be studied using the transfer matrix method introduced 
by Blote and Nightingale |^] , which takes advantage of the Fortuin-Kasteleyn represen- 
tation |^2| in terms of graphs of the partition function of the Potts model in order to 
reduce the dimension of the Hilbert space 0. In the Fortuin-Kasteleyn representation, 
the partition function (with no magnetic field) is 

-Z' = Tr + 

where uij = e^*^ — 1, is expanded as a sum over all the possible graphs Q (with s sites 
and 1{Q) loops) leading to the random cluster model: 

bij £ {0; 1} being the bond variables. Blote and Nightingale suggested to introduce a set 
of connectivity states which contain the information about which sites on a given row 
belong to the same cluster when they are interconnected through a part of the lattice 
previously built. A unique connectivity label rji = rj is attributed to all the sites i of 
such a cluster. In the connectivity space, \Z{m)) is a vector whose components are given 
by the partial partition function Z{m, {rii}m) of a strip of length m whose connectivity 
on the last row is given by {rii}m- The connectivity transfer matrix is then defined 
according to \Z{m + 1)) = Trn\Z{m)) and the partition function of a strip of length m 
becomes \Z{m)) = YVk=i T/c|Z(l)), where \Z{1)) is the statistics of uncorrelated spins. 



j j+u 

Figure 5: Transfer matrix along the disordered strip. 
The following physical quantities are measured: 

refined algorithm based on a loop representation of the partition function was proposed in Ref. |96| 
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i) The quenched free energy density is given (up to a ksT factor) by the Lyapunov 
exponent of the product of an infinite number of transfer matrices [97| 



(27) 



Ao(L) = hm — In 

m— >oo 777, 



(28) 



where | fo) is a unit initial vector. 

For a pure system, the central charge c is defined as the universal coefficient in the 
lowest-order correction to scaling of the free energy density /l of a strip of width 
L: 

fL = foo-^ + o(^), (29) 



6L2 



L4 



where the regular contribution is /oo = limi_»+oo fi- For a disordered system, c 
is defined in the same way from the finite-size behaviour of the quenched average 
free energy density /l, and numerically, since the strip widths available are small, 
we can only expect to measure effective central charges which depend on the dis- 
order strength, Ccs{g), and which would converge towards the true value c in the 
thermodynamic limit 96]. 



6L^ + "^^ ■ 



(30) 



ii) The spin-spin correlation functions in the time-direction (u) of the strip are cal- 
culated using an extension of the Hilbert space that allows to keep track of the 
connectivity with a given spin. For a specific disorder realisation, the spin-spin 
correlation function along the strip 



1 



q-l 



(31) 



is given by the probability that the spins along some row, at columns j and j + u, 
are in the same state and is expressed, in the absence of long-range order, in terms 
of a product of non-commuting transfer matrices: 



3"i + u 




(32) 



where | Aq) is the ground state eigenvector and T'^ is the transfer matrix in the 
extended Hilbert space. The operators gj and dj^^ realise the mapping between 
the two connectivity spaces. The correlations were computed on strips of varying 
widths and then averaged over many disorder realisations. 
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4 Analysis of numerical data in 2D 



4.1 Location of the random fixed point 

The transition line between ordered and disordered phases in the phase diagram starts at 
some point corresponding to the pure system and ends at another point where the critical 
properties are governed by the percolation universality class. Somewhere between, the 
random fixed point governs the critical behaviour of quenched randomness. Although 
this random fixed point is attractive, its precise location is an important preliminary 
step. Indeed, if the assumption of the existence of a unique stable random fixed point 
holds, one expects that the critical behaviour is asymptotically the same as the system 
is moved along the transition line. However, in finite systems, one generically has to 
deal with strong crossover effects due to the competition between the disordered fixed 
point and the pure and percolation fixed points, or to corrections to scaling linked to the 
appearance of irrelevant scaling variables. These latter effects are generally important 
in random systems and the corresponding corrections to scaling can be substantially 
reduced when one measures the critical exponents in the regime of the random fixed 
point, expected to be reached at the vicinity of the maximum of the effective central 
charge. 



0.9 
0.7 

0.5 



0.3 



G.. 
■ . ^ 

O.J 



10 



10 



G G- 



100 



L 



100 



1 ^ 

0.9^' 
0.8 



1000 



Figure 6: Finite-size scaling behaviour of the magnetisation (q = 
for different disorder amplitudes r (binary distribution), 
scaling are smaller close to r = 10 (taken from Picco |]28|| ' 



8, binary disorder) 
The corrections to 



Let us consider the finite-size behaviour of an observable Q measured at a deviation 
= K — Kc{g) from the critical point on some system of characteristic size L, in the 
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K 



Figure 7: Dependence of the effective central charge Ceg with the exchange coupHng K 
and the bond probabihty p for the 4-state Potts model. The maximum gives 
the location of the transition line and the absolute maximum corresponds to 
the optimal disorder strength. 

presence of disorder whose strength is measured by an amplitude g (ratio r between 
strong and weak interactions, probability p of non vanishing bond,. . . ). The variables t 
and play the role of relevant scaling fields (with positive RG eigenvalues yt = l/v 
and yi = 1 respectively), while close to the fixed point, disorder is supposed to be 
related to some irrelevant scaling variable with eigenvalue yg = —to < 0. At the fixed 
point there is no need that the irrelevant scaling field vanishes, so that one can write g* 
the corresponding disorder strength at the fixed point and the observable Q obeys the 
following homogeneity assumption in the scaling region 

Q(t, L-\g) = L--ViL'/% L-^ig - g*))- 

An expansion of the last variable (keeping the leading term only) along the critical line 
(i.e. varying g at Kdg)) gives 

Q(0, L-\g) = rQL--«(l + rg)(g - g*)L-^ + ...), 

where the F's are non-universal critical amplitudes. It is thus possible to fix g = g* in 
order to minimise the corrections to scaling (see figure 
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Figure 8: Phase diagram obtained from the condition of a maximum of the central charge 
for the dilute Potts model {q = 3, 4, and 8). The optimal disorder strength 
(here p*) is denoted by a larger circle and the full lines correspond to the 



single-bond effective medium approximation 100, 101, 102 1. 



From a finite-size scaling analysis (see section 4.3), the optimal disorder strength g 



is reached when a given quantity seems to be well fitted by a simple power-law (i.e. no 
bending in the the log-log plot). In the strip geometry, this value g* is found to coincide 



with the location of the maximum of the central charge along the critical line [43, 44 1, 
as a consequence of Zamolodchikov's c— theorem (see figure ^. 

In the literature, different types of disorder distributions have been considered. In 
the most studied case, when self-duality holds, the exact transition curve is known and 
the optimal disorder strength g* has to be found while moving one parameter only, which 
simplifies the task. If one defines the variables yij by e^''^ = -^^(e^'J — 1), the duality 
relation can be written |16] 



and self-duality is satisfied when the probability of each coupling equals the probability 
of its dual coupling, 

V{yi,)dy^,=V{y:^)dy*^, (33) 
that is if the distribution V{yij) is even. In the case of the bimodal distribution, the 
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self-duality point 

[eMKcir)) - 1] [exp(rKc(r)) - 1] = g , (34) 

corresponds to the critical point of the model if only one phase transition takes place in 
the system as rigorously shown in Ref. | |99| . 

In the case of dilution, self-duality does not work and the transition line has to be 
found numerically. The condition of a maximum of the central charge is again used as 
illustrated in figures and |8|. The result is in fair agreement with the effective medium 
approximation [|lO^, |lOl|, |TU^ . 

4.2 Temperature dependence 

According to their definition, the critical exponents can be obtained from a temperature- 
dependence study. Using for example the case of the order parameter in the low tem- 
perature phase, one can write 

M{t) = B^\tf{l + ...), t = Kc-K<0. 

1.0 
0.8 
^-0.6 
0.4 
0.2 
0.0 

0.0 0.1 0.2 0.3 0.4 0.5 

t 

Figure 9: Temperature dependence of the bulk and boundary magnetisation (g = 8, 
r = 10, binary disorder). 

The dots indicate the correction terms whose importance depend on the size of the 
system and on the distance from the transition temperature. Technically, one uses the 
definition of an effective temperature-dependent exponent, as illustrated in figures ^ and 

d in r t^o 

The precise value of the disorder strength g of course also infiuences the value of Pesit) 
(playing a role in the corrections to scaling as mentioned above) but asymptotically the 
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Figure 10: Temperature dependence of the exponents associated to the bulk and bound- 
ary magnetisation (q = 8, r = 10, binary disorder). 

limit t ^ should be independent on g, since there is only one fixed point governing 
the disordered system. 

4.3 Finite-size scaling 

One of the simplest method to extract critical exponents (once the critical temperature is 
known) is probably standard Finite-Size Scaling. On a finite system, the physical quan- 
tities cannot exhibit any singularity. They can be written as a singular term corrected 
by some scaling function which depends on the characteristic sizes of the problem, the 
correlation length ^ and the size of the system L. In the case of the order parameter den- 
sity for example we get mL{T) = \K — Kc\^ f{L/^,). The function f{x) of course depends 
on the geometry, but at the critical point Kc, the following behaviour is obtained: 

MUK,) ~ L-^/r (35) 

L—*oo 

Here, the ratio f3/v is precisely the magnetic scaling dimension Xg-. An example is shown 
in figure |ll] for the 8— state Potts model. From the slopes of the curves, the values 
-i/u = 1.686(17), 13/u = 0.152(4) and v = 1.005(30) can be obtained Q. The results 
here are interesting as reference values that we shall compare with more sophisticated 
techniques later. 

4.4 Short-time dynamics scaling 

It is commonly believed that universality can be found only in equilibrium stage of long- 
time regime in numerical simulations. For a magnetic system far from criticality, e.g. in 
the high temperature phase, suddenly quenched to the critical temperature, a universal 
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Figure 11: Finite-size scaling analysis of magnetisation, susceptibility and effective tran- 
sition temperature shift for the critical self-dual binary disordered 8— state 
Potts model. 



dynamic scaling behaviour emerges already within the short-time regime, according to 
a simple generalisation of the homogeneity assumption for the order parameter |54, 55 



56, 



M(t, r, L, Mo) = 6-^/^M(6i/'^t, 6"^, b'^L, b'^'^Mo), (36) 

where z is the dynamic exponent (dependent on the choice of algorithm), t = \K — Kc\ is 
the deviation from the critical point, Mq is the initial magnetisation with the associated 
scaling dimension xq, and r is the time (measured in MC sweeps). In the thermodynamic 
limit, L DO, and at criticality, t = 0, the expected evolution is given by M(r, Mq) = 
T^^f^^ f{MQT~^°^^) and allows the computation of the critical exponents. The main 
interest of short-time dynamics scaling is that it is not affected by critical slowing down, 
since only early time stages of the simulation are involved. 

4.5 Conformal mappings 

Monte Carlo simulations of two-dimensional spin systems are generally performed on 
systems of square shape while transfer matrix computations are done in strip geome- 
tries. In the following, we consider such a system of size 2N x L, and call u and v the 



corresponding directions (figure 12). The order parameter correlations between a point 
close to the surface, and a point in the bulk of the system should, in principle, lead 
to both surface and bulk critical exponents. Practically, it is not of great help for the 
accurate determination of critical exponents, since 
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i) strong surface effects (shape effects) occur which modify the large distance power- 
law behaviour, 

ii) the universal scaling function entering the correlation function is likely to display 
a crossover before its asymptotic regime is reached (system-dependent effect). 




Figure 12: Monte Carlo simulations of the 2d Ising model inside a square of 101 x 101 
lattice sites (10^ MCS/spin, Swendsen-Wang cluster algorithm). The figure 
shows the correlation function between a point close to the surface {wi = i) 
and all other points w in the square. The sketch on the right specifies the 
notations. 

One can proceed as follows: systems of increasing sizes are successively considered, 
and the correlations are computed along u- (parallel to a square edge considered as the 
free surface) and u-axis (perpendicular to this edge). The order parameter correlation 
function for example is supposed to obey a scaling form which reproduces the expected 
power-law behaviour in the thermodynamic limit: 




where Xa and are, respectively, the bulk and surface order parameter scaling dimen- 
sions. The scaling functions /^^ depends on the geometry, but have to satisfy asymptotic 
expansions including corrections to scaling e.g. /^^ (^) ~ 1 + const (^)^ + . . . in the 
boundary region u — > L. 

Equations ( p7| ) are not very useful for the determination of critical exponents, since 
at least five unknown quantities appear, and the correction terms in /^^ (x) may have a 
large amplitude, resulting from the significance of finite-size corrections. Nevertheless, 
if the power-law fit is limited to the linear regions on a log-log scale (almost one decade 
in figure p^), one gets correct estimations of the critical exponents, as long as condition 
ii) is fulfilled. 

At a critical point, scale invariance coupled with rotation and translation invariance 
also implies covariance under local scale transformations, i.e. conformal transforma- 
tions Q. The conformal transformations are those general coordinate transformations 
which preserve the angle between any two vectors. They leave the metric invariant up 
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Figure 13: Monte Carlo simulations in the case of the two-dimensional Ising model. 

a) Log-log plot of the order parameter correlation function perpendicu- 
lar to the surface (v-direction) and parallel to this surface (u-direction) 
for system sizes 101 x 101 (+, 5.10*^ MCS/spin) and 501 x 501 (o, 2.10^ 
MCS/spin). The correlations are computed between a point wi = i 
close to the surface and points w = iv, and w = u + i respec- 
tively. The corresponding estimations for + and 2x^ are writ- 
ten in the figure and the fit has been done in the range indicated by a 
dashed line. The figures correspond to power-law fits for the size 501^. 

b) Conformal rescaling of the perpendicular correlations for the two sizes 
(see below). The rescaled correlation function exhibits a true power-law in 
the whole range of variation of the new variable, and the value for Xa- + x^. 
is improved. The solid lines correspond to the theoretical asymptotic form 
while the dotted line is the exact expression in the continuum limit. 



to a scale change. The conformal group includes the Euclidean group as a subgroup, 
the dilatations and the special conformal transformations. In any dimension these con- 
formal transformations are also called global conformal transformations because they 
map the infinite space onto itself. Consider a lattice model described in the continuum 
limit by a local theory defined by some action S. Conformal symmetry occurs when a 
local theory is scale invariant. In lattice systems, nearest-neighbour interactions ensure 
that physics is local and scale invariance holds at the critical point where a continuous 
limit description is allowed. For any local field (the energy density or the magnetisation 
for a spin system for example) the usual homogeneity assumption under a homogeneous 
rescaling R — > 6R 

is extended to local transformations with a position dependent rescaling factor. In 
two dimensions conformal transformations are realized by the analytic functions in the 
complex plane (the conformal group is thus infinite-dimensional): z — > w{z) and equa- 
tion (38) is thus generalised to a covariance law of transformation of the (A^— point) 
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Figure 14: Conformal mapping of the infinite plane z = x+iy inside a finite square NxN 
with free edges or a strip of width L with periodic boundary conditions. 



correlators under conformal mappings: 



(39) 



Here, zi and Z2 are two points in the original complex plane and wi and W2 are 
the corresponding points in the transformed complex geometry under the mapping 
w{z). This transformation law is very helpful in numerical analysis, since simula- 
tions or numerical computations are always performed on finite systems of particular 
shape, depending on the technique used. The critical properties of an infinite system 
{4i{zi)(j){z2)) ~ \zi — Z2\~'^^'f' can thus be obtained by fitting the numerical data to 
the transformed conformal expression which usually deviates significantly from a simple 
power law (although the algebraic decay at criticality must be recovered asymptotically 
in the limit of an infinite system of course). The situation is schematically sketched in 



figure 14 

Among all the possible mappings, we shall here specify a few cases of interest: 
i) Mapping onto a cylinder: the logarithmic transformation 

w[z) = — mz = u + IV 
2tt 



(40) 



is well known to map the infinite plane onto a strip of finite width L with peri- 
odic boundary conditions and infinite length. This is the limit — > oo in the 
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rectangular geometry described above. Using the algebraic decay of the two-point 
correlator in the z— plane, one gets on the strip 



{(f>{0,0)4>iu,v)) 



T 



2 cosh 



/2ttu 



2 cos 



/'2Trv 



(41) 



In the long direction of the strip, at large distances it becomes an exponential 
decay 



{cP{0,0)cP{u,0))pkc= [^j exp(^ 



For sufficiently large strip widths, the transverse direction can also give some inter- 
esting results. Using the mapping w{z) = ^Inz, the half-infinite plane is mapped 
onto a strip with open boundaries in the transverse direction. If the boundary 
conditions are fixed (for example using an ordering surface field coupled to the 
order parameter) on one edge and free on the opposite edge (this is the meaning 
of the notation +f below), the transverse profile of the order parameter density is 
given by the conformal expression 



{cr{v))^f = const x 



L f irv 
- sm — 
vr \ L 



F 



cos 



TTV 

2L 



(43) 



The shape of the scaling function F{x) is asymptotically constrained by simple 



scaling, F{x) 
parameter). 



~ X 



(here, x^ is the boundary scaling dimension of the order 




Figure 15: Transverse profiles of the order parameter in an infinitely long strip with 
fixed-free boundary conditions. 

ii) Mapping onto a square: the Schwarz-Christoffel transformation 

w{z) = —F{z,k), z = sn^-^j (44) 
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maps the half-infinite plane z = x + iy{0<y< oo) inside a square w = u + iv of 
size N X N {—N/2 < u < N/2, < v < N) with free boundary conditions along 
the four edges. Here, F{z, k) is the elliptic integral of the first kind, sn {2Kw/N) 
the Jacobian elliptic sine, K = K{k) the complete elliptic integral of the first kind, 
and the modulus k is solution of K{k)/K{\/1 — k"^) = i. 

In the semi-infinite geometry, the two-point correlator is fixed up to an unknown 
scaling function (apart from some asymptotic limits implied by scaling). Fixing 
one point zi close to the free surface [zi = i) of the half-infinite plane, and leaving 
the second point exploring the rest of the geometry, Z2 = z, the following behaviour 
is expected: 

(0(zi),/.(z))i^~2/--''VM, (45) 

where the dependence on u) = of the universal scaling function •i/' is con- 

strained by the special conformal transformation and its asymptotic behaviour, 
^/^(u;) ~ Lo^'t', in the limit y ^ 1, is implied by scaling. 

Using the mapping (^^, the local rescaling factor in equation (^) is obtained, 
w'{z) = ^[(1 — •2^)(1 — k'^z^)]"^/'^ , and inside the square the two-point correlation 
function becomes (see Ref. fl^ ) 

{<i,{w{)(t>{w)),^. ~ {^[z\. (|1 - z\\l - k''z^\y^'^^ "%(c^), (46) 
^ . ' 

k{w) 

with z{w) given by equation (^^. This expression is correct up to a constant 
amplitude determined by k{wi) which is kept fixed, but the function il){uj) is still 
varying with the location of the second point, w. 

In order to cancel the role of the unknown scaling function, it is more convenient to 
work with a density profile in the presence of ordering surface fields. This is a one- 
point correlator whose functional shape in the half-infinite geometry is determined 
by scaling apart from some amplitude: 

(^(■2))ioo = constxy-^- (47) 

and it maps onto 

(cj(u;))sq. = constx [k(w)]~'''" (48) 
where the function k{w) again comes from the mapping. 

The case of the mapping of the infinite complex plane inside a square with periodic 
boundary conditions was for example used in Ref. [^] . 

Other mappings can be convenient. Our choice here was motivated by the fact that 
Monte Carlo simulations are usually performed on samples of square shapes. On the 
other hand, the strip (with free or fixed boundary conditions) is the natural geometry 
generated in transfer matrix calculations. 
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Figure 16: Profiles of the order parameter in a square with fixed boundary conditions. 
4.6 Impact of rare events and non self-averaging 

As we mentioned, the physical properties have to be averaged over many samples pro- 
duced with a given probability distribution. One can typically encounter two opposite 
situations, depending on the quantities under interest. The average order parameter for 
example is defined as a sum of quantities affected by randomness, 

i 

while the correlation function 

^ {^0\ Sj { n n)d,+JAo) 

k=j 

essentially depends on a product of non-commuting matrices whose elements are de- 
termined by the disorder distribution. These two types of quantities definitely exhibit 
different properties as functions of the number of samples used to sample the probability 
distribution. When computing mean values, it is clear that the accuracy of the results 
for a given number of disorder realizations, does not have the same behaviour in the case 
of sums or of products of random variables. Consider as an example the sum and the 
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product of random variables Aj taken from a binary distribution, 

n n 

1=1 1=1 

and compute the moments (rescaled by a power 1/p) [{T,x)p]^^^ and averaged 
over N realizations of the Aj's (here we choose for this example n = 50, so there are 
some 10^^ configurations). Numerical results are given in the table |^ for some choice of 
parameters. 









Sa 








Ha 




TV 




p = 5 


p = 20 


p = 50 


p=l 


p — 5 


p = 20 


p = 50 


10^ 


0.990 


1.060 


1.228 


1.305 


1.054 


1.529 


2.042 


2.188 


10^ 


1.014 


1.049 


1.156 


1.298 


1.060 


1.292 


2.143 


2.630 


10^ 


1.010 


1.047 


1.161 


1.327 


1.054 


1.303 


2.480 


3.467 


10^ 


1.010 


1.047 


1.160 


1.322 


1.055 


1.301 


2.511 


3.857 



Table 3: Comparison between the moments of a sum and of a product of random vari- 
ables distributed according to a bimodal probability distribution, as a function 
of the number of realizations. 

It is particularly clear that the values obtained from the sum S converge rapidly (the 
variations between results obtained at different number of realizations correspond to a 
statistical noise). This is due to the fact that E is normally distributed, while in the 
case of the product IT, we note a continuous increase of the numerical estimate of a given 
moment as the number of samples increases and this effect is especially pronounced for 
high moment order p. This is due to the fact that the distribution of the II's is log-normal 
and thus there exist some events which have a dominant role in the average, but which 
are so rare that they are not scanned by a poor statistics which would essentially explore 
the region of typical values. In order to be more precise in the distinction between typical 
and average value, we rewrite In 11 as a sum of random variables. In IIa = J2i 1^ which, 
according to the central limit theorem, has a Gaussian distribution in the limit of large 
number of draws. The typical value corresponds to the maximum of the probability 
distribution, Iltyp = e^'^'^ where the Gaussian is centered, and clearly differs from the 
average value 11 = e^"'^. 

The same situation occurs when computing the spin-spin correlation function 

{a{0)a{u)),, = G;^. (49) 

Since it is almost log- normal (see figure 0), the logarithm of Ga{u) is self-averaging 
and the average \n.G(j{u) as well as higher order moments are well behaved. A cumulant 
expansion is thus convenient to reconstruct the average Ga{u) through 

G„{U) = el"! («)+ I (Iri^ (M)-ln Ga {uf) + 
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Figure 17: Top: Probability distribution of the spin-spin correlation function (8-state 
Potts model with bimodal disorder). The insert shows V{Ga{u)) (with a 
very long tail on the right), while it is shown on a logarithmic scale in the 
main frame, where one can notice the shape which is close to a log-normal 
distribution. Bottom: Reconstruction of the correlation function from the 
moments of its logarithm. 



This is a test (see figure 17) which proves that the probability distribution is sufficiently 
well scanned with the large numbers of realizations used in this run. 
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5 Numerical results and comparison with perturbative 
expansions in 2D 

5.1 Regime q > 4 

5.1.1 Randomness induces a second-order regime 

Although there is no perturbation result for g > 4, we shall start in the regime of first- 
order transition of the pure model. A first step was to prove that even for large numbers 
of states per spin, the transition was rounded to become continuous. That was done 
by different authors Chen et al studied the free energy barrier AF{L), defined 

from the energy histogram V{E) in Monte Carlo simulations according to e"^^^*-^^ = 

-fmax/ -Pwelh 

with Pmax given by the maximum of V{E) and P„eii corresponding to the 
value at the bottom of the well separating the two coexisting phases. They showed 
that the energy barrier AF{L) = —2ao,d.L'^~^ vanishes in the thermodynamic limit 
(figure ^)where iTo.d. is the order-disorder interface tension between the two possibly 
coexisting phases. 

The dynamics of the Monte Carlo simulations leads to compatible conclusions. The 
energy autocorrelation time te is indeed exponentially large (with the system size) when 
a non vanishing order-disorder interface tension (To.d. exists, 

while it is only increasing as a power law at second-order transitions, 

Te ~ L^, 

with a dynamical exponent z which strongly depends on the algorithm used. This latter 



mea- 



situation is indeed observed |Q in the disordered 8— state Potts model (figure 18). 

Cardy and Jacobsen on the other hand used transfer matrix calculations [41| 
suring the free energy density fi whose corrections to scaling behave at first-order tran- 
sitions like /l ~ /oo + ©(L^'^e-^/^). Plotting then A(L) = ln(/L - f^o) + dlnL vs the 
strip width L should give asymptotically a straight line with a slope given by the in- 
verse correlation length 1/^. In the presence of randomness, the curve corresponding 
to the 8— state Potts model indicates a diverging correlation length as expected at a 
second-order phase transition (figure p^). 



5.1.2 Comparison between finite-size scaling and conformal mappings 

Conformal mappings provide quite efficient techniques for the determination of critical 
exponents. The validity of such an approach is nevertheless restricted to systems where 
scale invariance, translation invariance and rotation invariance do hold. This require- 
ment is not obviously fulfilled in random systems, since disorder breaks the symmetries. 
Hopefully, one may expect that after averaging over many disorder realizations, one is 
led to some effective system for which these symmetries are restored. This assumption 
can be checked from numerical simulations. Studying the critical behaviour using an 
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Figure 18: Top: Evolution of the free energy barrier with the size of the system for 
the 8— state Potts model with binary disorder (at small disorder strength). 
The dotted line in the inset shows the case of the pure system (taken from 
Chen, Ferrenberg and Landau Ref. [^]). The dotted line is a guide for the 
eyes. Bottom: Power-law behaviour of the energy autocorrelation time in the 
8— state Potts model with binary disorder. 



independent technique, namely finite-size-scaling which can safely be supposed to give 
the right results, we then compare to the exponents deduced from various mappings 
which are assumed to work. The comparison was performed carefully in the case of the 
8— state Potts model with binary disorder, and the technique was then applied in the 
regime q > A, and even to asymptotically large g's The FSS results, which are 

considered here as the reference results, are shown in figure In strip geometries, 
the correlation functions in the long direction and the order parameter profile in the 
transverse direction (figure 20) lead to compatible results. 

In the square geometry, the fit of the correlation function must be done along some 
curves inside the square where the scaling variable w remains constant, which implies 
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Figure 19: Evolution of the corrections to scaling to the free energy in a strip geometry 
(taken from Cardy and Jacobsen |4l[|). 



that the function ^l){u)) also remains constant and a simple power-law fit is then needed, 

{a{wi)a{w))^^ ~ A^\k{w)\-''\ 



where stands for the amplitude which contains {wi is fixed and chosen equal to 
imaginary unit wi = i). A difficulty occurs due to discretization of the lattice Estimation 
of the correlation function along the continuous curve oo = const, is done using a Taylor 
expansion from the data taken at 4 neighbour points on each plaquette. This explains 



the non monotonic behaviour of the error bars in figure 21. 

The case of the density profile is better to achieve, since it involves no unknown 
scaling function, but a single constant amplitude (it is a one-point correlator). All the 
points inside the square enter the power-law fit 



(cr(u;))gq ~ const x \n{w)\ 

and make it more accurate (and there is no need of any expansion around the lattice 
points). The number of points {N'^) being so large, we can even forget uncertainties on 
each point and simply take as error bar on the resulting exponent the standard deviation. 

The results for the magnetic scaling dimension measured using Finite-Size Scaling 
techniques, compared to the mappings onto strips or square geometries, are compared 
in table ^ in the case of the 8— state Potts model with a self-dual binary probability 
distribution of coupling strengths. First we note that the transition is second-order, 
even in the regime g > 4, as predicted by the Imry-Wortis agreement. More important 
for the following is the fact that the agreement between different techniques is quite fair 
and leads to the conclusion that conformal techniques can be applied here, in spite of 
the lack of the symmetry properties which should in principle be required. This is due 
to the fact that we are interested in average quantities. The system thus becomes, on 
average, translationnaly, rotationally invariant, as well as scale invariant at the critical 
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q=S, binary disorder 
0.1 1 

sin nv/L 

Figure 20: Fixed-free BC order parameter critical profile for the self-dual binary disor- 
dered 8— state Potts model. The behaviour close to the fixed surface gives 
access to the bulk scaling dimension (the corresponding curve is shown in dot- 
ted line and the conformal expression is shown in full line while the symbols 
correspond to strips of various widths). 



point. This is an important point, because the use of conformal mappings is more 
accurate than standard FSS methods, and the comparison between different schemes in 
the 2 < q < 4 regime will require a great accuracy, as already noticed in table ^. 

5.2 Regime q < 4 

5.2.1 Tests of replica symmetry 

First of all, the fitting procedure has to be validated. From the exponential decay of the 
average spin-spin correlation function, the exponent x'^ is deduced and presented as a 



function of q for the case of a binary disorder in Fig. 22. These results were first reported 
by Cardy and Jacobsen in Ref. The agreement with the third order expansion in 
equation (|20| ) is extremely good especially in the region where the expansion is supposed 
to be valid when q is not too far from the Ising model value q = 2. The quality of 
the data confirms the reliability of the averaging procedure (table |5|). Even close to the 
marginally irrelevant case of the Ising model where logarithmic corrections are known to 
be present for some quantities, we note that the numerical data are quite satisfactorily in 
agreement with the perturbative results. The agreement is made better by the absence 
of logarithmic corrections for the average correlation function at q = 2, and we will see 
that this observation is no longer true in the following study of other moments. 

The question of a possible breaking of replica symmetry in disordered systems is very 



controversial and far from being settled, especially in spin glasses (see e.g Refs. |103, 104 , 



105, 106, 107|). In the context of disordered Potts ferromagnets, the question was first 



addressed by Dotsenko et al |25|. In order to test between Replica Symmetry and Replica 



Symmetry Breaking schemes, Dotsenko et al performed a second order expansion of the 
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Figure 21: Fitting the correlation function (top) and the order parameter profile (bot- 
tom) in a square geometry. The different curves (top) correspond to different 
values of the variable u. This is the case g = 8, with a binary distribution, 
at the optimal disorder strength estimated by the maximum of the central 
charge in strip geometry. 



exponent of the second moment of the spin-spin correlation function decay in both cases 
(Equations (|2^)). MC simulations were first performed at q = 3 but were not completely 
conclusive, although in favour of Replica Symmetry: the perturbation expansion leads 
to x'^2(3) = 0.1176 and x'^2(3) = 0.1201 according to equations (p^, while previous 
numerical results lead to 0.113(1) [||], 0.1140(5) ||, 0.116(1) || and 0.119(2) ||. 

Conclusive results for different values of q were then obtained using transfer matrices. 
Close to g = 2, the proximity of the marginally irrelevant Ising FP will surely alter the 
data, as a reminiscent effect of the logarithmic corrections present exactly at g = 2 for 
the second moment []T7| , |108| . Too large values of q on the other hand are not very helpful 
in order to check perturbation expansions which break down when one explores higher 
values of the expansion parameter. One thus has to balance between these two extreme 
situations and the comparison between numerical data and perturbation results should 
be conclusive around g = 3 or slightly below. The TM technique thus appears to be 
well adapted, since it is capable to deal with non integer values of q. The comparison is 
shown in Fig. ^ for the bimodal probability distribution and the results are also given 
in table |^ In the convenient domain for the test, around (7 = 3, results are written 
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Figure 22: Top: Scaling dimension of the order parameter (binary disorder) compared 
to the third order expansion of Dotsenko and co-workers |^0|. The scahng 
dimension corresponding to the pure model is shown for comparison. Bot- 
tom: Exponent of the second moment of the spin-spin correlation function 
as a function of the number of states of the disordered Potts model (binary 
disorder). The comparison is done with Replica Symmetry and Replica Sym- 
metry Breaking scenarios [^]. The agreement with the RS result is quite 
good around g = 3. When q is close to 2, the discrepancy can be attributed 
to the weak relevance of disorder. We indeed used a simple exponential fit as 
can be expected at a stable disordered FP, but at q = 2, one knows from Lud- 
wig's results that logarithmic corrections must be added. These corrections 
can also influence the vicinity of g = 2 in a numerical approach. 
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for the 8— state Potts model with binary disorder 



Technique 



Quantity 



Scaling dimension Ref. 



Standard techniques 

t— dependence Mb{t) 

FSS MbjKc) 

FSS (g(0)g (L/2)) 

Short-time dynamics Mi,(t) 



Conformal mappings 

Periodic strip 
Free BC square 
Fixed-free strip 
Fixed BC square 



0.151(1) 
0.153(1) 
0.159(3) 
0.151(3) 



|33| 
|28, 
|34| 

l"l 



(a(0)a(u;)),q 



0.1505(3) 
0.152(3) 
0.150(1) 
0.1503(1) 




Table 4: Comparison between temperature-dependence, FSS and short-time dynamics 
scaling results for the magnetic exponent and the scaling dimension x'^ 

at the random fixed point deduced from the logarithmic and the Schwarz- 
Christoffel mappings for the 8— state Potts. In the first line, the exponent 
P' deduced from the temperature dependence is close to x'^, since the value of 
the correlation length exponent is found very close to 1 (i^' ~ 1.01(1)). We note 



also that in all these references but Ref. |34], a binary distribution of disorder 
was used. 



in bold face. The agreement with Replica Symmetry is quite convincing. The results 
for the exponent associated to the average energy-density correlations confirm this 
conclusion. 

5.2.2 Multiscaling 

The multiscaling behaviour of the spin-spin correlation functions is noticeable in the 
dependent set of exponents of the reduced moments 

In Ref. |46U , exhaustive computation of 50 different moments in the range < p < 5 

were performed in the strip geometry, and the associated scaling dimensions followed 

from a semi-log fit \n {a{0)a{u))P vs Inn, followed by an extrapolation to L — > oo. The 

numerical results were compared to the first order expansion of Ludwig and to the second 

order expansion in the RS scheme in Eq. (^). The second order result is clearly very 

good up to values of p close to 3 and then breaks down as already noticed by Lewis . 

An alternate presentation of the results (used e.g. by Ludwig |p^) is given b y the 

scaling dimension of the moment of the correlation function itself, {a{0)a{R))P (not 
i/p 

the reduced function {a{0)a{R))P ). The scaling dimension is thus simply px'^p{q), 
hereafter denoted by X'^p(q). An example, with g = 3, is shown in Fig. ^ where we 
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Expansion (bol) TM result (Rcf. M) 



2 

2.25 
2.5 
2.75 
3. 

3.25 
3.5 



0.12500 
0.12800 
0.13051 
0.13269 
0.13465 
0.13653 
0.13845 



0.1252(3) 

0.1282(5) 

0.1307(7) 

0.1328(9) 

0.1347(11) 

0.1364(13) 

0.1379(14) 



Table 5: Comparison of the numerical results for the magnetic scaling dimension (bi- 
modal probability distribution) x'^ with the third order expansion of Dotsenko 
and co-workers ||2^. The error bars systematically contain the analytical value. 

have shown the results obtained with the bimodal and the continuous self-dual prob- 
ability distributions at the optimal disorder amplitude as well as the dilution case at 
optimal dilution. Once again, we find a fair agreement between the numerical data and 
the perturbative result which confirms universality, i.e. the exponent associated to a 
given moment of the correlation function does not depend on the detailed probability 
distribution of the coupling strengths. 

5.2.3 Probability distribution of correlation functions 

In Ref . [jl9| , Ludwig presented a remarquable discussion of the spin-spin correlation func- 
tion probability distribution. The relevant information on the large distance behaviour 
is encoded in the multiscaling function H{a), which is simply the Legendre transform 
of the set of independent scaling indexes X'^p(q). Setting dX'^p{q) = adp, this function 
is simply obtained by 

H{a)=X'^p{q)-ap. 

The geometrical interpretation of this Legendre transform follows from the relation ^ = 

—p where a is defined by ^'^Qp^'^^ = ol. The scaling dimension x'^p{q) is obtained on 
the plot of H{a) by the intercept of the tangent of slope —p with the abscissa axis. 
An example of multiscaling function H{a) deduced from the numerical data with the 
bimodal probability distribution is shown in Fig. 

In this section, we follow Ludwig's arguments and report a numerical study of the 
correlation function probability distribution in the cylinder geometry. 

According to the results of the previous section, the moments of the spin-spin corre- 
lation function along the strip asymptotically behaves as follows: 

GHu) = {a{id)a{u))P ~ Bpe'^Kv (51) 

and are defined in terms of the probability distribution V[G{u)]: 

GHuj= [\g{u)V[G{u)]GP{u). (52) 
Jo 
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q 

2.25 
2.5 
2.75 
3. 

3.25 
3.5 



Perturbative results TM result 



0.12213 

0.12002 

0.11854 

0.11761 

0.11718 

0.11723 



0.12229 
0.12067 
0.11997 
0.12011 
0.12110 
0.12304 



(Ref. |46||) 



0.1204(5) 
0.1194(8) 
0.1185(10) 
0.1177(12) 
0.1172(14) 
0.1169(16) 



(Rcf- M) 



2.5 

2.75 

2.5 



1.006 
1.013 
1.023 



1.000 
1.000 
1.000 



i.oo(iy 

1.01(1) 
1.02(1) 



Table 6: Decay exponent of the second moment of the spin-spin correlation function 
compared to Replica Symmetry and Replica Symmetry Breaking expressions of 
The results written in bold face correspond to the range of values 

The second part of the 



Eqs. (H). 

of q where the agreement is particularly satisfactory. 



table presents Jacobsen's results |47] for the exponent of the average energy 
correlation function. 



Following Ludwig, we introduce the variable Y{u) = —lnG{u) and write G^{u) 
^-pY{u)^ Using the identity V[G{u)]dG = V[Y{u)]dLY and equations (|l|) and (| 
one obtains 

/•no ^ 

dy(n)P[y(^x)]e-P^(") ~ Bpe-^^'- 



L 



which leads to the expression of the probability distribution by inverting the Laplace 
transform (5 > 0): 

V[Y(u)] = lim / dp5„e ^ \- 5^pJ . 

2m 5^0 Js-ioo 

The amplitude Bp is assumed to be smoothly dependent on p (this can be checked 
numerically), and its dependence can be forgotten with respect to the exponential, since 





average 


typical 






TM 


x',„ TM 


2.5 


1.006 


1.00(1) 


1.023 1.02(1) 


2.75 


1.013 


1.01(1) 


1.051 1.04(2) 


3. 


1.023 


1.02(1) 


1.090 1.06(3) 



Table 7: Decay exponent of the average and typical energy-energy correlation functions: 
comparison between perturbative expansion from Ref. |27] and transfer matrix 
computation (from Jacobsen |47]). 
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Figure 23: Comparison of the multi-fractal exponents (moment of the correlation func- 
tion {a{0)a{u))P) with the second order expansion of Lewis in the RS 
scheme for both the bimodal and the continuous probability distribu- 
tions. 



it only introduces a correction when 2ttu/L — > oo. Let us define the function h{p) = 
X'^p — 2^|t/L ^- III the large distance limit 2ttu/L — > oo, the integral can be evaluated by 
the saddle-point approximation at the minimum po of h{p): 



^ X' 



pa 



Yju) 
2ttu/L 



Instead of Y{u), we define the scaled variable a = ' saddle point value at pQ 

only depends on this variable /i(po) = H{a). We thus obtain the probability distribution 



V[Y{u)] ~ exp 
or, using V[Y{u)]dY = V{a)da, 

Via) 



27ru 



jj ( Yj^) \ 
L \2ttu/lJ 



2ttu 



L 



exp 



27rn 



H(a) 



(53) 



(54) 



The multi-fractal function contains the essential information on the probability distri- 
bution. A correction to the leading behaviour given by the saddle-point approximation is 
needed here to improve the data collapse onto a single multiscaling function. If we expand 
the function h{p) close topo) h{p) ~ H{a) + \h" [pojip—pof' , with h"{po) > we obtain, 
instead of Eq. (|53|), the following result for the probability distribution P[y(ti)] [|109(| : 



V[Y{u) 



2'nu 

IT 



"1/2 



exp 



2iru 



( Yju) 
L " \2ttu/L 



H 



(55) 
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and a correction appears in V{a): 

Via) 



/27rn\i/2 



27rn 



H{a) 



(56) 



which enables to extract H{a) at fixed a by fitting the probabihty distribution to the 
expression 



lnV{a) 



1 , 2ttu 
- In 

2 L 



const 



2ttu 

~17 



H{a). 



(57) 



It is shown in Fig. ^ where the probabihty distribution of the spin-spin correlation 
function was obtained after collecting the results over 96 000 disorder realizations in 50 



classes All the data collapse onto a single multiscaling function H{a), for different 
distances u and this is even the case with still good accuracy for different strip widths 



or different probability distributions |46 
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6 Conclusion and summary of the main results 



The two-dimensional Potts model is the ideal framework to test the influence of quenched 
randomness on phase transitions and critical behaviour. It exhibits a second-order phase 
transition completely characterised by conformal invariance when the number of states 
per spin is lower or equal to 4 and a first-order transition above. The transition line is 
exactly known, and it is easy to build in the random case probability distributions of 
coupling strengths which preserve the self-duality relation. 

With respect to these advantages, many results concerning the eff^cct of a weak dis- 
order were obtained during the last decade using perturbations expansions around the 
pure fixed point. Different solutions were considered, first replica symmetric solutions 
where the symmetry between all the replicas is supposed to be preserved in the renor- 
malization equations, then the case of a spontaneous breaking of the replica symmetry 
was also studied perturbatively. 

Numerical studies were also performed from different sides. First of all, Monte Carlo 
simulations coupled to finite-size scaling analysis, then transfer matrices and sophisti- 
cated graph and loop algorithms coupled to extensive use of conformal mappings in 
order to extract the values of the critical exponents with a pretty good accuracy. All 
the results were in a fair agreement with the perturbative expansions close to the Ising 
model limit and concluded in favour of replica symmetric scenarios. The multiscaling 
properties of the order-parameter and energy density were also analysed, and a charac- 
terisation of the probability distributions of the correlation functions was made in terms 
of universal multiscaling functions. The regime q > 4 where the pure model exhibits 
a first order transition was also extensively studied, but did not display any particular 
features compared to the regime q < 4 in the presence quenched randomness. 

The most important question which remains open up to now is the identification 
of the conformal field theories which could describe the random fixed point. This is a 
delicate program, since in the replica limit of coupled models, the theory is not unitary 
(with a central charge which evolves continuously with the value of q). We may hope 
some progress in this direction for the near future, which would undoubtedly achieve a 
considerable progress in the understanding of two-dimensional disordered systems. 
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